Nighttime Lights Analysis

Author

Data Lab

Published

April 18, 2024

Nighttime lights have become a commonly used resource to estimate changes in local economic activity. This section shows where nighttime lights are concentrated across Syria and observed changes over time.

Data

We use nighttime lights data from the VIIRS Black Marble dataset. Raw nighttime lights data requires correction due to cloud cover and stray light, such as lunar light. The Black Marble dataset applies advanced algorithms to correct raw nighttime light values and calibrate data so that trends in lights over time can be meaningfully analyzed. From VIIRS Black Marble, we use data from January 2012 through present—where data is available at a 500-meter resolution.

Map of nighttime lights

We first show a map of nighttime lights. Most of the country is dark, with lights concentrated within cities.

Code
## Load boundaries
adm0_sf <- read_sf(file.path(boundaries_dir, 
                             "ner_adm_ignn_20230720_ab_shp",
                             "NER_admbnda_adm0_IGNN_20230720.shp"))

## Load/prep raster
prep_r <- function(year_i){
  r <- rast(file.path(ntl_dir, "individual_rasters", "annually",
                      paste0("VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",year_i,".tif")))
  r <- r %>% mask(adm0_sf)
  r[][r[] == 0] <- NA
  r[] <- log(r[] + 1)
  r[] <- log(r[] + 1)
  return(r)
}

r_2012 <- prep_r(2012)
r_2013 <- prep_r(2013)
r_2014 <- prep_r(2014)
r_2015 <- prep_r(2015)
r_2016 <- prep_r(2016)
r_2017 <- prep_r(2017)
r_2018 <- prep_r(2018)
r_2019 <- prep_r(2019)
r_2020 <- prep_r(2020)
r_2021 <- prep_r(2021)
r_2022 <- prep_r(2022)
r_2023 <- prep_r(2023)

## Make map
pal <- colorNumeric(c("yellow", "orange", "red"), unique(c(r_2012[],
                                                           r_2013[],
                                                           r_2014[],
                                                           r_2015[],
                                                           r_2016[],
                                                           r_2017[],
                                                           r_2018[],
                                                           r_2019[],
                                                           r_2020[],
                                                           r_2021[],
                                                           r_2022[],
                                                           r_2023[])),
                    na.color = "transparent")

leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addRasterImage(r_2012, colors = pal, opacity = 1, group = "2012") %>%
  addRasterImage(r_2013, colors = pal, opacity = 1, group = "2013") %>%
  addRasterImage(r_2014, colors = pal, opacity = 1, group = "2014") %>%
  addRasterImage(r_2015, colors = pal, opacity = 1, group = "2015") %>%
  addRasterImage(r_2016, colors = pal, opacity = 1, group = "2016") %>%
  addRasterImage(r_2017, colors = pal, opacity = 1, group = "2017") %>%
  addRasterImage(r_2018, colors = pal, opacity = 1, group = "2018") %>%
  addRasterImage(r_2019, colors = pal, opacity = 1, group = "2019") %>%
  addRasterImage(r_2020, colors = pal, opacity = 1, group = "2020") %>%
  addRasterImage(r_2021, colors = pal, opacity = 1, group = "2021") %>%
  addRasterImage(r_2022, colors = pal, opacity = 1, group = "2022") %>%
  addRasterImage(r_2023, colors = pal, opacity = 1, group = "2023") %>%
  addLayersControl(
    baseGroups = paste0(2012:2023),
    options = layersControlOptions(collapsed=FALSE)
  )

Maps of change in nighttime lights

Pixel level

Code
r_load_month_pc <- function(month_i){
  r_b1 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                         paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_01",".tif")))
  r_b2 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                         paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_02",".tif")))
  r_b3 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                         paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_03",".tif")))
  r_b4 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                         paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_04",".tif")))
  r_b5 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                         paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_05",".tif")))
  r_b6 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                         paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_06",".tif")))
  r_e <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                        paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t",month_i,".tif")))
  
  r_b <- r_b1
  r_b[] <- (r_b1[] + r_b2[] + r_b3[] + r_b4[] + r_b5[] + r_b6[]) / 6
  
  r_e[] <- (r_e[] - r_b[]) / r_b[] * 100
  r_e[][r_b[] <= 0.5] <- NA
  
  r_e[][r_e[] >= 100] <- 100
  r_e[][r_e[] <= -100] <- -100
  
  r_e <- r_e %>% crop(adm0_sf) %>% mask(adm0_sf)
  
  return(r_e)
}

r_23_07 <- r_load_month_pc("2023_07")
r_23_08 <- r_load_month_pc("2023_08")
r_23_09 <- r_load_month_pc("2023_09")
r_23_10 <- r_load_month_pc("2023_10")
r_23_11 <- r_load_month_pc("2023_11")
r_23_12 <- r_load_month_pc("2023_12")
r_24_01 <- r_load_month_pc("2024_01")
r_24_02 <- r_load_month_pc("2024_02")
r_24_03 <- r_load_month_pc("2024_03")

## Make map
r_values <- unique(c(r_23_07[],
                     r_23_08[],
                     r_23_09[],
                     r_23_10[],
                     r_23_11[],
                     r_23_12[],
                     r_24_01[],
                     r_24_02[],
                     r_24_03[]))

pal <- colorNumeric(c("red", "white", "green"), r_values,
                    na.color = "transparent")

leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addRasterImage(r_23_07, colors = pal, opacity = 1, group = "2023-07") %>%
  addRasterImage(r_23_08, colors = pal, opacity = 1, group = "2023-08") %>%
  addRasterImage(r_23_09, colors = pal, opacity = 1, group = "2023-09") %>%
  addRasterImage(r_23_10, colors = pal, opacity = 1, group = "2023-10") %>%
  addRasterImage(r_23_11, colors = pal, opacity = 1, group = "2023-11") %>%
  addRasterImage(r_23_12, colors = pal, opacity = 1, group = "2023-12") %>%
  addRasterImage(r_24_01, colors = pal, opacity = 1, group = "2024-01") %>%
  addRasterImage(r_24_02, colors = pal, opacity = 1, group = "2024-02") %>%
  addRasterImage(r_24_03, colors = pal, opacity = 1, group = "2024-03") %>%
  addLayersControl(
    baseGroups = c("2023-07",
                   "2023-08",
                   "2023-09",
                   "2023-10",
                   "2023-11",
                   "2023-12",
                   "2024-01",
                   "2024-02",
                   "2024-03"),
    options = layersControlOptions(collapsed=FALSE)
  ) %>%
  addLegend("bottomright", pal = pal, values = r_values,
            title = "% Change in NTL<br>relative to<br>Jan-June 2023",
            opacity = 1,
            labels = c("< -100", "-50", "0", "50", "> 100")
  )

ADM level

Code
adm3_sf <- read_sf(file.path(boundaries_dir, 
                             "ner_adm_ignn_20230720_ab_shp",
                             "NER_admbnda_adm3_IGNN_20230720.shp"))

adm3_sf <- adm3_sf %>%
  select(-date)

adm3_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm3_monthly.Rds"))

adm3_df <- adm3_df %>%
  mutate(baseline = (date %in% c(ymd("2023-01-01"),
                                 ymd("2023-02-01"),
                                 ymd("2023-03-01"),
                                 ymd("2023-04-01"),
                                 ymd("2023-05-01"),
                                 ymd("2023-06-01")))) %>%
  group_by(ADM3_PCODE) %>%
  mutate(ntl_sum_baseline = mean(ntl_sum[baseline %in% 1])) %>%
  ungroup() %>%
  
  mutate(ntl_pc = (ntl_sum - ntl_sum_baseline) / ntl_sum_baseline * 100,
         ntl_change = ntl_sum - ntl_sum_baseline) %>%
  select(date, ADM3_PCODE, ntl_pc, ntl_change, ntl_sum, ntl_sum_baseline) %>%
  
  filter(date >= ymd("2023-01-01")) 

adm3_df$ntl_sum_baseline_dates <- adm3_df$ntl_sum
adm3_df$ntl_sum_baseline_dates[adm3_df$date >= ymd("2023-07-01")] <- NA

make_ntl_date <- function(adm3_df, ym_i){
  adm3_df[[paste0("ntl_sum_", ym_i)]] <- adm3_df$ntl_sum
  adm3_df[[paste0("ntl_sum_", ym_i)]][adm3_df$date != ymd(paste0(ym_i, "-01"))] <- NA
  return(adm3_df)
}

adm3_df <- adm3_df %>%
  make_ntl_date("2023_07") %>%
  make_ntl_date("2023_08") %>%
  make_ntl_date("2023_09") %>%
  make_ntl_date("2023_10") %>%
  make_ntl_date("2023_11") %>%
  make_ntl_date("2023_12") %>%
  make_ntl_date("2024_01") %>%
  make_ntl_date("2024_02") %>%
  make_ntl_date("2024_03")

adm3_df$ntl_pc[adm3_df$ntl_pc >= 100]  <- 100
adm3_df$ntl_pc[adm3_df$ntl_pc <= -100] <- -100

## Add sparkline
ntl_date_color <- "black"

adm3_df_line <- adm3_df %>%
  filter(date >= ymd("2023-01-01")) %>%
  arrange(date) %>%
  split(.$ADM3_PCODE) %>% 
  map_df(~{
    l_ntl_all <- sparkline(.x$ntl_sum,
                           type='bar',
                           barColor="orange",
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = 'orange', 
                           highlightSpotColor = 'orange')
    l_ntl_base <- sparkline(.x$ntl_sum_baseline_dates,
                            type='bar',
                            barColor="blue",
                            chartRangeMin = 0,
                            chartRangeMax = 8,
                            width = 100,
                            height = 50,
                            tooltipChartTitle = "Nighttime Lights",
                            highlightLineColor = 'blue', 
                            highlightSpotColor = 'blue')
    l_ntl_2023_07 <- sparkline(.x$ntl_sum_2023_07,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2023_08 <- sparkline(.x$ntl_sum_2023_08,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2023_09 <- sparkline(.x$ntl_sum_2023_09,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2023_10 <- sparkline(.x$ntl_sum_2023_10,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2023_11 <- sparkline(.x$ntl_sum_2023_11,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2023_12 <- sparkline(.x$ntl_sum_2023_12,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2024_01 <- sparkline(.x$ntl_sum_2024_01,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl_2024_02 <- sparkline(.x$ntl_sum_2024_02,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
      l_ntl_2024_03 <- sparkline(.x$ntl_sum_2024_03,
                           type='bar',
                           barColor=ntl_date_color,
                           chartRangeMin = 0,
                           chartRangeMax = 8,
                           width = 100,
                           height = 50,
                           tooltipChartTitle = "Nighttime Lights",
                           highlightLineColor = ntl_date_color, 
                           highlightSpotColor = ntl_date_color)
    l_ntl <- spk_composite(l_ntl_all, 
                           l_ntl_base)
    data.frame(l_ntl = as.character(htmltools::as.tags(l_ntl)))
  }, .id = 'ADM3_PCODE') 

## Merge data with sf object
adm3_df <- adm3_df[!is.na(adm3_df$ADM3_PCODE),]
adm3_sf <- adm3_sf %>%
  right_join(adm3_df, by = c("ADM3_PCODE")) %>%
  left_join(adm3_df_line, by = "ADM3_PCODE")

adm3_sf <- adm3_sf %>%
  dplyr::mutate(popup = paste0("<h4>", ADM3_FR, "</h4>",
                        "<b>NTL ", substring(date, 1, 7), ":</b> ", round(ntl_sum, 2), "<br>",
                        "<b>NTL Baseline:</b> ", round(ntl_sum_baseline, 2), "<br>",
                        "<b>Trends in NTL since Jan 2023</b><br><br>",
                        "<span style='color: blue;'>Blue</span> = Jan - June 2023<br>",
                        l_ntl))

## Make map
r_values <- unique(adm3_sf$ntl_pc)

pal <- colorNumeric(c("red", "white", "green"), r_values,
                    na.color = "transparent")

line_color <- "black"
fill_opacity <- 0.8
opacity <- 0.5

add_deps <- function(dtbl, name, pkg = name) {
  tagList(
    dtbl,
    htmlwidgets::getDependency(name, pkg)
  )
}

leaflet() %>%
  addProviderTiles(providers$OpenStreetMap.HOT) %>%
  addPolygons(data = adm3_sf[adm3_sf$date %in% ymd("2023-07-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2023-07") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2023-08-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2023-08") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2023-09-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2023-09") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2023-10-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2023-10") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2023-11-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2023-11") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2023-12-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2023-12") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2024-01-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2024-01") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2024-02-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2024-02") %>%
  addPolygons(data = adm3_sf[adm3_sf$date == ymd("2024-03-01"),], 
              fillColor = ~pal(ntl_pc), 
              color = line_color, opacity = opacity, fillOpacity = fill_opacity, 
              label = ~lapply(popup, HTML),
              group = "2024-03") %>%
  onRender("function(el,x) {
      this.on('tooltipopen', function() {HTMLWidgets.staticRender();})
    }") %>%
  addLayersControl(
    baseGroups = c("2023-07",
                   "2023-08",
                   "2023-09",
                   "2023-10",
                   "2023-11",
                   "2023-12",
                   "2024-01",
                   "2024-02",
                   "2024-03"),
    options = layersControlOptions(collapsed=FALSE)
  ) %>%
  addLegend("bottomright", pal = pal, values = r_values,
            title = "% Change in NTL<br>relative to<br>Jan-June 2023",
            opacity = 1,
            labels = c("< -100", "-50", "0", "50", "> 100")
  ) %>%
  add_deps("sparkline") %>%
  browsable()

Percent change in nighttime lights from 2023 to 2024

Below we show the percent change in nighttime lights for each month in 2024 relative to the same month in 2023.

Code
## Load data
adm0_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm0_monthly.Rds"))
adm1_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm1_monthly.Rds"))

## Add dummpy ADM1 variable so function will more easily work on both datasets
adm0_df$ADM1_FR <- ""

## Function to prep data as percent changes
prep_month_pc <- function(adm_df){
  
  adm_long_df <- adm_df %>%
    filter(date >= ymd("2023-01-01"),
           date <= ymd("2024-06-01")) %>%
    select(ADM1_FR, date, ntl_sum) %>%
    mutate(month = date %>% month,
           year = date %>% year) %>%
    pivot_wider(id_cols = c(ADM1_FR, month),
                values_from = ntl_sum,
                names_from = year) %>%
    filter(!is.na(`2024`)) %>%
    mutate(pc = (`2024` - `2023`) / `2023` * 100)
  
  return(adm_long_df)
}

adm0_pc_df <- prep_month_pc(adm0_df)
adm1_pc_df <- prep_month_pc(adm1_df)

## Maps
adm0_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  theme_classic2() +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across Niger from 2023 to 2024")

Code
adm1_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  theme_classic2() +
  theme_custom +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across regions from 2023 to 2024") +
  facet_wrap(~ADM1_FR,
             scales = "free_y")